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ABSTRACT 
The  problem  of  maximum  range  atmospheric  reentry  for  an  orbiting 
lifting  glider  was  treated  by  Bryson  and  Denham  by  the  method  of 
"steepest  ascent"  <,    The  same  problem  is  undertaken  here  by  a  method 
of  differential  corrections  developed  by  Faulkner ,  This  method  makes 
use  of  a  Newton-Raphson  type  iteration  based  on  paths  which  satisfy 
the  Euler-Lagrange  equations,  A  comparison  of  results  is  madeD  show= 
ing  large  differences  in  control  variable  history „  and  longer  range 
for  the  path  obtained  by  differential  corrections .  The  problem  was 
characterized  by  a  sharp  " ridge '!  in  the  domain  of  the  starting  values 
of  the  adjoint  variables  and  the  effect  of  this  on  the  convergence  of 
both  methods  is  discussed.  Finally ,  the  difficulty  of  choosing  initial 
approximations  for  the  starting  values  of  the  adjoint  variables  is 
discussed,  and  a  method  is  presented  for  obtaining  these  from  a  nominal 
path  as  the  first  step  in  the  computer  routine. 


1 
A,  E0  Bryson  and  W,  F.  Denham ,  A  steepest-ascent  method  for  solv= 
ing  optimum  programming  problems  (Ratheon  report  BR»1303),  Raytheon 
Company,  Bedford,  Massc,  1961. 
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INTRODUCTION 

A  problem  of  optimum  control  may  be  described  very  roughly  as 
the  determination  of  the  decision -making  or  control  function  which 
will  result  in  the  largest  possible  value  of  the  final  "payoff"  from 
some  continuing  process „ 

The  process  might  be  illustrated  as  shown  in  figure  1<,19  where 
P  is  the  "payoff" B  something  whose  value  depends  on  the  terminal  point 
or  the  path  history ,  and  Q  is  a  set  of  given  constraints  which  must  be 
satisfied „ 


Given  initial  point. 


t  =  T 

P  =  maximum 

Q  =  0 


Figc  1 e1   An  Optimum  Path 


The  methods  of  solution  commonly  used  now  were  foreseen  at  least 
forty  years  ago,  but  the  volume  of  computation  required  limited  their 
application  to  only  the  simplest  of  problems „  In  the  past  ten  years 
there  has  been  an  explosion  of  interest  due  to  the  availability  of 
high-speed  digital  computers  on  the  one  hand,,  and  on  the  other  hand 
to  the  urgent  need  for  solutions  in  such  applications  as  space 
trajectories  o  The  "state  of  the  art"  is  such  that  no  method  is  com- 


pletely  general ,  and  a  great  diversity  of  approaches  is  employed 9 
reflecting  the  absence  of  a  firm  framework  of  mathematical  theory 
concerning  nonlinear  differential  equations D 

In  this  thesis  a  comparison  is  made  between  two  methods  for  cal- 
culating optimum  angle  of  attack  programs  for  a  lifting  unpowered 
vehicle  to  attain  maximum  range „  starting  from  circular  orbital 
velocity  at  a  height  of  3°0  thousand  feet  above  the  surface  of  the 
earth.  The  problem  was  chosen  because  it  is  typical  of  a  class  of 
problems  in  which  there  is  now  great  interest 9  and  because  two  of 
the  foremost  workers  in  the  field  of  optimizations,  Arthur  E„  Bryson 
and  W,  Fo  Denham„  had  published  a  clear  and  well-documented  solution 
by  the  method  of  "steepest  ascent" „ 

The  comparison  was  suggested  by  Stanley  Ross,  while  one  of  the 
writers  was  engaged  in  summer  field  work  under  the  sponsorship  of 
the  Lockheed  Missile  and  Space  Company,  John  V,  Breakwell  and  George 
Leitmann  generously  gave  many  hours  of  their  time  to  general  dis- 
cussions and  made  specific  suggestions  which  were  very  helpful » 
Professor  Frank  D.  Faulkner  of  the  U8  S0  Naval  Postgraduate  School 
not  only  laid  out  the  general  form  of  the  problem  and  overcame  the 
difficulties  that  arose ,  but  also  the  method  of  solution  used  was 
the  one  developed  by  him  from  the  basic  method  of  differentials  of 
Bliss, 


Chapter  I 
GENERAL  FORMULATION  OF  A  PROBLEM  OF  OPTIMUM  CONTROL 

1.1  The  equations  of  motion,  the  control  variables  and  the  con- 
straints . 

It  will  be  helpful  to  define  an  optimum  control  problem  in 
somewhat  more  exact  terms  than  those  given  in  the  introductory 
remarks . 

Suppose  we  are  given  a  system  whose  behavior  is  described 
by  a  system  of  differential  equations  which  we  shall  call  the 
equations  of  motion.  We  will  suppose  that  wherever  any  of  these 
differential  equations  is  of  higher  than  first  order,  we  have  de- 
fined such  additional  variables  as  necessary  so  as  to  obtain  n 
first  order  differential  equations  of  the  form 

S^    =    *i^l»    ®2*     •••»    Sji9    p-i9    Pop     •••«    Pm»     W 

1      "~      J-  £>       lC|        c   o   a   ji      n 

where  the  variables  s.  are  the  state  variables ,  the  p^  are  the 
control  variables,  t  is  the  independent  variables  and  the  dot 
denotes  differentiation  with  respect  to  the  independent  variable, 
It  is  assumed  that  the  f.  are  of  class  C". 

It  will  be  convenient  for  what  is  to  follow,  if  matrix 
notation  is  adopted  at  this  point.  The  equations  of  motion  may 
be  written  as  the  single  matrix  equation 
(1.1)  S  ■  F(S,  P,  t) 


where 


S  = 


n 


,  an  nxl  matrix  of  the  dependent  state 
variables 


F  ■ 


n 


,  an  nxl  matrix  of  known  functions  of  the 
control  variables,  the  state  variables , 
and  possibly  the  independent  variable 


P  = 


m 


,  an  mxl  matrix  of  control  variables , 
which  we  are  free  to  choose  as  functions 
of  t. 


In  addition,  there  are  k  constraints  or  end  conditions,  of  the  form 
e.(S,t)  =  0,  which  may  define  values  of  the  state  variables  at  the 
end  point  t  =  T,  or  may  be  of  integral  form  such  as 


Kz  h 


dt  =  constant, 


In  the  latter  case,  we  shall  define  an  additional  variable  sn+-j_ 
so  that 


s     =  f. 
n+1    i 


with  the  initial  and  end  conditions 


sn+l  <°>  =  ° 


'n+1  (T)  =  the  given  constant. 


These  may  be  written  as 


(1.2)     E  = 


,  a  kxl  matrix  of  constraint  functions 


We  wish  to  discover  the  particular  control  variable 
program  which  will  drive  this  system  from  a  given  initial  point 
to  some  terminal  point  which,  of  the  set  of  all  such  points  which 
satisfy  the  terminal  constraints,  gives  a  maximum  value  to  some 
function,  M(S,t),  of  the  state  variables;  that  is,  for  the  ter- 
minal point  t  =  T,  possibly  unknown, 

M(S,t)T  =  M(T)  =  max. 

1.2  Bounded  control  variables  and  inequality  constraints. 

The  permissible  domain  of  the  control  variables  may  likewise 
be  limited.  Some  examples  of  restrictions  on  the  control  for  com- 
mon physical  systems  might  be: 


f  pi  ^  P  max, 


f   P  dt^I  max, 


P  =  Pmax  for  0<t<tx  and 

P  =  0     for  t1<t<T. 
Bounded  control  variables  appear  frequently  in  real  physical 
systems.  For  solving  these  problems  it  may  be  useful  to  re- 
place the  inequalities  by  mathematically  equivalent  equalities „ 
This  is  accomplished  by  the  introduction  of  suitable  new  vari- 
ables in  the  manner  due  to  Valentine  (ref.l).  For  instance, 
for  a  constraint 

0£p£P  max 
one  may  define  the  real  variable  g,  so  that 

P(pmax  -  P)  -  g2  =  0. 
This  procedure  transforms  the  inequality  to  an  equality,  and  we 
adjoin  this  equation  to  the  set  (l.l).  The  same  procedure  is 
applied  in  the  case  of  inequality  constraints  or  end  conditions. 

We  will  generally  require  that  all  state  variables  be  speci- 
fied initially,  but  this  is  not  necessary.  One  or  more  of  these 
may  be  unspecified,  that  is  "free". 

The  general  form  of  these  problems,  then,  is  the  two-point 
boundary  value  problem.  With  the  stipulation  that  the  functional 
M  is  to  be  maximized  (or  minimized)  they  become,  depending  on  the 
formulation,  the  classical  problem  of  Bolza  or  problem  of  Mayer 
of  the  calculus  of  variations.  Since  to  maximize  some  quantity, 
say  u,  is  equivalent  to  minimizing  the  quantity  minus  u,  we  will 
henceforth  consider  that  the  term  "maximum"  includes  both  cases 


unless  otherwise  specified. 

1.3  Degenerate  problems. 

Some  important  observations  may  be  made  here.  First,  as 
pointed  out  by  Faulkner  (ref.2)  and  Breakwell  (ref.3)»  the  case 
in  which  all  of  the  functions  f .  involve  any  function  of  a  par- 
ticular control  variable  only  linearly  are  said  to  be  degenerate . 
These  degenerate  problems  typically  are  of  the  "bang-bang"  type 
in  which  the  control  changes  dis continuously.  The  neglect  of 
induced  drag  in  aerodynamic  problems  for  instance  often  results 
in  bang-bang  control. 

Second,  the  constraints  must  not  involve  the  control  variables, 
or  the  problem  is  overspecified.  If  the  control  program  is  to  be 
such  as  to  produce  an  extreme  of  the  functional  M,  then  it  may  not 
simultaneously  satisfy  any  prescribed  constraint  at  any  isolated 
point  of  the  path,  else  the  solution  is  generally  discontinuous. 

1.4  The  adjoint  variables. 

We  shall  introduce  a  set  of  Lagrange  multipliers,  as  yet  un- 
determined 

U=  fl'  U2 \| 

a  lxn  matrix  of  adjoint  variables.  The  term  'adjoint"  is  used, 

since  they  will  be  chosen  to  be  solutions  to  the  system  of  differ- 
ential equations  which  is  adjoint  to  equations  for  the  variations 
of  the  given  system. 

Suppose  we  write  (1.1)  as 


S  -  F  =  0 
and  integrate  its  product  with  the  adjoint  vector  between  the 
limits  of  the  independent  variable.  For  convenience,  we  will 
call  these  limits  zero  and  T,  and  we  have, 

/  U(S  -  F)  dt  =  0. 
o 

The  corresponding  variational  equation  is 

/TU(^S  -  F  6S  -  F  f?)   dt  =  0. 
o        s     P 

Now  this  is  integrated  by  parts  to  eliminate  from  the  integrand 
the  variations  of  the  state  variables. 

*k+ 

k- 


(1.3)  fi^l   T  =/TR  U+UF  )5S  +U   F6Pldt-jllu   f1      k    dt 


where  Fs  and  F_  are  the  matrices  of  partial  derivatives, 


(l.M 


F.  = 


af1/ps1  .  .  .  2>f1/asn 


9fn/3s^  •    •    •  ^^*n'^sn 


,  an  nxn  matrix, 


FP  = 


af1/ap1  .  .  .  a^/apj^ 


afn/3pi  •   .  .  3*jf2Pm 


,  an  nxm  matrix, 


and  t  is  a  symbol  for  any  point  or  set  of  points  where  S  or  U 
is  discontinuous. 

Now,  if  U  is  chosen  as  a  solution  to  the  equation 


8 


(1.5)  U  +  U  F  =  0, 

5 


then  (3)  becomes 

—  —  T       T  +  + 

(1.6)      )u6s]   =  /  (U  F  S?)   dt  -  L  Fuf]  k  dt 


This  sequence  of  operations  defines  the  system  which  is  adjoint 

to  the  variational  equations  of  (1.1),   OS  -  F  6S  =  F  6?   , 

s      p 

and  forms  the  1  x  n  matrix,  or  vector 


...       . 

U  =   u.,  u0  .  .  .  u 
J_  1  2       n 


Hereafter,  no  distinction  will  be  drawn  between  a  vector  in 
k-space  and  alxk  or  kxl  matrix,  or  a  corresponding 
column  or  row  of  a  matrix. 
1.5  Green's  formula  and  the  Euler-Lagrange  equations. 

Equation  (1.6)is  the  fundamental  formula  relating  the  varia- 
tions of  the  end  values  of  the  state  variables  with  the  variations 
of  the  control  variables.  Note  that  the  formula  shows  the  ter- 
minal value  of  the  adjoint  variable  to  be  the  sensitivity  of  the 
terminal  value  of  the  corresponding  state  variable  to  a  change  in 
control.  This  interpretation  will  be  discussed  further  in  con- 
nection with  the  transversal  conditions  (sect.  1.  9  ). 

Equation  (1.6)  is  often  called  Green's  formula  (one  form). 
See  Coddington  and  Le Vinson  (ref.  ^)  page  86.  The  last  term  on 


the  right  applies  only  at  points  of  discontinuity  of  F,  and 

o 

since  these  are  discontinuities  in  S,  they  are  ■corners". 
Hereafter,  it  will  be  assumed  that  the  solution  does  not  have 
a  corner  unless  it  is  otherwise  mentioned.  In  this  case  the 
last  term  drops  out.  The  fundamental  lemma  of  the  calculus  of 
variations  states  that,  given  a  control  program  P(t)  such  that 
the  constraints  E  are  satisfied,  then  if  the  coefficient  of  6P 
in  (1.6)  does  not  vanish  for  some  solution  U  to  the  adjoint 
it  is  possible  to  construct  a  variation  such  that  the  constraints 
are  still  satisfied,  but  the  function  M  exceeds  the  value  on  the 
first  path.  A  proof  of  this  lemma  is  included  in  refs,  (6)  and  (10) 
Thus  for  M  to  have  at  least  a  stationary  value  M* , 
(1.7)  UF  =  0.     (The  Suler  equation) 

Equation  (1.5)  defines  the  system  of  equations  which  is  adjoint 
to  the  variations  of  the  equations  of  motion,  (1.1).  The  adjoint 
equations  and  equation  (1.7)  will  be  called  the  Euler-Lagrange 
equations.  It  will  be  shown  that  these,  together  with  the  con- 
straints and  the  equations  resulting  from  the  transversal  condi- 
tions comprise  a  complete  set  which  determines  the  optimum  solution. 

1.6  The  maximum  principle  of  Pontryagin. 

The  functions  f.  are  often  composed  of  terms  which  are  func- 
tions of  the  state  variables  alone  and  terms  which  contain  also 
the  control  variables.  Suppose  we  collect  terms  of  the  first  type 


10 


into  an  n  x  1  matrix  we  will  call  V,  and  those  of  the  second  type 
into  an  n  x  1  matrix  G.  Then  (1.1)  becomes 

S  =  V  +  G. 
Now  consider  these  matrices  to  be  vectors  in  an  n-dimensional 
hyperspace  of  the  state  variables,  and  a  particularly  fruitful 
geometric  interpretation  of  (1.7)  appears.  Let  us  consider  equa- 
tion (1.7) 

U  •  t   =  0 

P 

as  the  condition  that 

U  •  G  =  max 
P 

that  is,  let  us  require  that  the  vector  G  be  everywhere  so  dir- 
ected that  the  projection  of  G  upon  U  is  maximized.  The  control 
program  which  maintains  this  relationship  (while  satisfying  the 
constraints)  is  optimum  at  least  locally.  This  statement  of  the 
condition  is  due  to  the  Russian  mathematician  Pontryagin,  and  it 
leads  to  particularly  enlightening  geometric  formulations  in  many 

■hi 

problems.  Note  that  the  shape  of  the  domain  of  G  then  gives  in- 
sight into  the  behavior  of  the  control  program. 

In  figure  1.2  a  hypothetical  domain  of  G  is  shown  by  the 
dashed  line.  The  scalar  quantity  U  •  G  is  maximized  when  G  has 
a  maximum  projection  on  U  at  all  times.  Figure  1.3  shows  the 
situation  where  more  than  one  G  vector  satisfies  the  maximum 
condition.  This  generally  yields  a  discontinuity  of  the  control, 
hence  a  "corner". 
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The  functions  gi  will  be  called  the  "forcing  functions". 
It  is  the  vector  whose  components  are  these  functions  which  must 
be  directed  so  as  to  maximize  its  projection  on  the  adjoint  vec- 
tor. The  equation  or  equations  which  express  this  condition  in 
the  form 

Pi  =  h(U,  S,  P,  t) 

are  obtained  by  ordinary  calculus  from  (1.7).  They  are  sometimes 
called  "steering  equations". 

1.7  The  Hamiltonian. 

From  the  foregoing,  one  sees  that  the  product  UF  has  special 
significance  when  U  is  chosen  according  to  (1.5) •  For  convenience, 
we  shall  define  this  quantity  as  the  Hamiltonian  function,  H. 

H  =  UF. 

The  maximum  principle  becomes 


(1.3)  H  =  max  , 

P 


and  the  general  condition  for  an  extreme  of  H  is 

(1.9)  H  =  0, 

P 

which  is  equivalent  to  (1.7). 

For  any  one  extremal,  if  t  does  not  appear  explicitly  in  the 
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equations  of  motion,  the  Hamiltonian  is  constant ,     This  may  be 
shown  readily,   for 

H  ■  UF  +  UF 

=  UF  +  UF  S  +  UF  P  +  UFX 
s  p  t 

=  (U  +  UFjF  +  UTP  +  UF.  . 

S  p  "L 

But  U  +  UF    =  0  by  equation  (1.5;,  and  UF     =0  by  (1.7),  hence 
s  p 

(1.10)         H  =  UFt  =  Ht 

In  a  great  many  problems ,  the  independent  variable  does  not 
appear  explicitly  in  the  equations  of  motion,  and  the  Hamiltonian 
becomes  a  constant  of  the  system.  In  the  next  section  it  will  be 
shown  that  in  this  case  if  the  final  value  of  the  independent  vari- 
able, T,  is  free  (that  is,  it  is  not  the  quantity  to  be  optimized 
and  is  not  constrained)  then  H  =  0  everywhere  on  a  given  path, 
provided  the  Euler-Lagrange  equations  (1.5)  and  (l.b)  are  every- 
where satisfied  on  the  path. 

1.8  Systems  which  are  linear  in  the  state  variables. 

Consider  the  functional 

T      T 
j[uF  dt  =f0E   dt  =<p(T)  . 

If  the  differential  equations  of  motion  are  linear  in  the  state 
variables,  then  the  functional  does  not  involve  the  state  variables 
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and  an  absolute  maximum  may  be  obtained.  If  the  differential 
equations  are  not  linear,  then  there  are  no  corresponding  gen- 
eral proofs,  ana  we  can  look  only  for  a  stationary  value  and  a 
local  extreme  for  the  functional. 

A  system  of  differential  equations  which  is  linear  in  the 
state  variables  with  constant  coefficients  may  be  written  as 

S  =  A  S  +  Q 
where  S  is  the  nxl  matrix  of  equation  (1.1),  Q  is  an  nxl  matrix 
whose  ith  component  is  the  forcing  function  of  the  control  vari- 
able in  the  ith  direction,  and  A  is  an  nyji    matrix  of  constants. 
Now,  in  analogy  with  the  derivation  of  equation  (1.6), we  may  form 
the  integral 

J    U(S  -  AS  -  Q)dt. 
o 

Integrating  by  parts 

_      T  T 

|USj      =     /  p  +   UA  )S  +  UqI  dt. 
oo 

Taking  U  to  be  the  solution  to  the  adjoint  equation  U  +UA  =  0,  we 

have  the  form  of  the  Green-Lagrange  formula 

T  T 

(1.11)  [US]     =    /UQdt. 

o       o 

Suppose  that  we  have  by  some  means  discovered  a  set  of 
initial  values  for  the  adjoint  variables  and  integrated  forward 
until,  at  some  time  T, 
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H  1.  the  curve  satisfies  the  constraints  on 

the  state  variables  and  the  control  vari- 
able, 
H  2.  the  forcing  function  Q  is  such  that  UQ  is 

maximum  everywhere  on  the  curve,  and 
H  3»  the  terminal  value  of  the  kth  component 
of  the  adjoint  vector  is  unity,  and  all 
other  components  are  zero  at  t  =  T, 
In  this  case,  Faulkner  (ref.  5)  shows  that  the  curve  furnishes  an 
absolute  maximum  for  the  kth  state  variable  under  the  given  initial 
values  and  the  constraints.  Akheizer  (ref,  6,  sect.  2.15)  and 
Edelbaum  (ref.  7»  sects.  1.1  and  1.2)  discuss  conditions  for 
strong  and  weak  extremes  in  this  and  in  nonlinear  cases. 

1.9  The  transversal  condition. 

The  extremals  have  been  discussed,  as  the  solution  curves 
which  satisfy  the  system 

(1.12) 


c 

0 

S  - 

F  =  0 

1  "  + 

UF     = 
s 

0 

H     » 
\     P 

=  UF     = 
P 

■    0 

Consider  the  system  already  discussed,  and  the  corresponding 
Green-Lagrange  variational  equation  (1.6), 


T  T  t  +  T 

Iu6sl       =   /(UF  6P)dt  -  £  [uf1     k  dt,    =f  6h  dt. 

LJoo         P  kLJt-k° 

k 
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We  require  that  the  tern  on  the  left,  evaluated  at  the  starting 
point  be  zero.  If  s.(o)  is  specified,  then  6s.  (0)  =  0.  If  s^ 
is  not  specified,  then  we  shall  choose  the  particular  solution 
to  the  adjoint  so  that  U.(0)  =0.  It  is  necessary  that  the  last 
term  on  the  right  be  zero,  else  changing  dt  will  offer  a  better 
trajectory.  There  are  variations  for  which  b?   is  not  "small"  and 

rT 

J   SH  dt  <0, 
o 

but  for  first  order  effects,  the  only  condition  on  bs   for  an  ex- 
tremal is  that 

(u  6s)T  =  0. 
The  form  of  the  equation  does  not  depend  on  the  constraints »   hence 
we  will  define  as  admissible  paths  those  which  satisfy  the  constraints, 
Since  the  constraints  must  be  satisfied  at  t  =  T,  the  endpoint  S(T) 
must  lie  on  the  n  -  k  dimensional  surface 

E(S,T)  =  0 
Hence,  for  an  admissible  differential  at  the  endpoint, 


de.  ■ 

i 


Z.A   ae./as.ds.  +  3e./dtdT    =0 


or  in  matrix  notation, 
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(1.13) 


dE  ■  (EsdS  +  EtdT)T  =  0 


where  dE  and  E^  are  of  dimension  kxl,  E  is  kxn  and  dS  is  nxl. 
Es  may  be  considered  to  be  the  gradient  of  E  in  n-space.  If 
Es  is  bordered  on  the  right  by  the  column  E+,  this  may  similar- 
ly be  considered  to  be  the  gradient  of  E  in  the  n+1  space  (of 
the  state  variables  and  t )  which  we  will  call  VE  . 


VE*  = 


de-^/ds^   .  .  .ae-|_/asn  ae-j^/at 


•  V  9 


Now  if  dS*  is  taken  to  be 


dS*  = 


ds. 


ds, 

dT 


n 


(1.13)  may  be  rewritten  as 


dE  =  (vE*dS*)T  =  0. 


In  the  same  manner,  we  may  have  the  gradient  of  M  in  n+1  space, 


VM*  = 


9K/as1  3K/5s2  .  0  .  dM/ds 


n  ^/atj 


and  the  augmented  adjoint  vector 

U* 


■  £  u2  •  •  •  \  (-h3 
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The  vector  U*  must  lie  in  the  manifold  spanned  by  the  k+1  vectors 
Vek,  VM  ,  hence  it  may  be  expressed 


U*  =  C 


VE* 

VM* 


=  C  Z* 


where  C  is  a  lx  (k+1)  row  vector  of  constants  and  Z*  is  the 
(k+l)  x  (n+l)  matrix  obtained  by  bordering  VE*  below  withVM*. 
Since  M  is  assumed  to  be  independent,   Z*  has  rank  k+1.  The 
transversal  condition  may  be  stated  as  the  condition  that  the 
matrices  VE*,  Z*  and  Y* 


Y*  = 


Z* 
U* 


,  a  (k+2)  x  (n+l)  matrix 


all  have  rank  k+1 . 

As  a  simple  example,  suppose  we  have  the  state  variables 
w,  x,  y  and  z,  with  all  initial  conditions  specified  and  for 
t=*T   w  constrained,  x  to  be  maximized,  y,  T  and  z  free.  We  may 
write 


w 


X 


10     0     0     0  constraints 

0     10     0     0  maximum 

u,    Un    Uo    u^   -H  adjoint  variables 
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a  3x5  matrix  which  must  be  of  rank  2  at  t  =  T.  Hence 

u  =  u^  =  -H  =  0  at  t  =  T. 

In  the  development  above,  free  use  was  made  of  material  from 

11 

a  research  paper  of  Faulkner.   When  each  e.  and  M  involves  only 

one  variable  of  the  set  (S,t),  the  condition  may  be  summarized  as 
f  olloitfs : 


(1.14) 


s  (T)  free  implies  u.(T)  =  0 


T  free 


implies  H(T)  =  0 
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Chapter  II 

SOLUTION  OF  A  GENERAL  OPTIMUM  PROGRAMMING  PROBLEM  USING 

DIFFERENTIAL  CORRECTIONS 

2.1  The  missing  constants  of  the  adjoint  set. 

It  has  been  shown  that  if  relations  (1,12)  and  (1.14)  are  satis, 
fied,  then  the  resulting  path  generally  furnishes  an  optimum  to  the 
state  S(T).  Finding  this  particular  path  depends  upon  finding  a 
particular  set  of  constants,  which  are  the  initial  values  of  the 
unknown  adjoint  variables.  This  being  true,  suppose  for  the  moment 
that  the  means  exist  for  correcting  a  set  which  is  "pretty  good", 
until  the  final  "perfect"  set  is  obtained.  That  is,  we  visualize 
in  the  domain  of  U(0),the  adjoint  vector  evaluated  at  t=0,  the 
point  U*(0)  whose  coordinates  are  these  "perfect"  starting  values. 
If  U(0)  =  U*,  then  integration  of  the  first  two  of  equations  (1.11) 
with  the  control  chosen  according  to  the  third  of  equations  (1.12) 
will  produce  the  optimum  path.  At  some  t=T  on  this  path,  the  con- 
straints will  be  satisfied,  the  desired  variable  will  be  maximiz- 
ed, and  simultaneously  (1.12)  and  (1.14)  will  be  satisfied. 

In  Chapter  III  a  method  will  be  given  that  will  allow  us  to 
correct  from  an  arbitrary  point  U'(0)  toward  U*(0)  provided  U'(0) 
is  within  some  unknown  region  R  in  the  vicinity  of  U*(0).   The 
difficulty  is  to  find  some  point  U(0)  which  lies  within  R.  The  so- 
called  Direct  Method  of  solution  used  here  reduces  the  problem  to 
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the  solution  of  a  two-point  boundary-value  problem  involving  the 
extremals.  Consequently,  the  missing  boundary  conditions  must 
somehow  be  guessed.  In  general,  the  guess  must  be  within  some 
region  R,  if  the  corrective  program  is  to  converge.  In  the  Grad- 
ient, or  Steepest  Ascent  methods  of  Kelley  and  Bryson  a  nominal 
trajectory  is  guessed  and  the  two-point  boundary-value  problem 
is  solved  by  correcting  along  paths  xrtiich  are  not  extremals.  This 
method  too  has  its  disadvantages.  For  some  problems  the  optimum 
path  may  never  be  reached,  although  the  correction  program  "con- 
verges", that  is,  satisfies  the  convergence  criterion.  For  the 
problem  treated  in  this  thesis,  a  path  was  obtained  by  the  direct 
methods  given  here  which  attained  a  final  value  of  the  variables 
to  be  maximized  well  beyond  that  of  the  "optimum"  path  obtained 
by  the  method  of  Steepest  Ascent. 

2.2  Nominal  paths. 

An  advantage  of  the  Gradient  Method  is  the  following.  In 
physical  optimization  problems  we  often  have  some  idea  beforehand 
what  kind  of  control  program  is  likely  to  come  fairly  close  to  the 
optimum  one.  The  only  "guessing"  involved  in  the  Gradient  Method 
is  to  choose  such  a  program  and  construct  a  nominal  path  which 
connects  the  given  initial  constraints  with  the  given  terminal 
constraints.  Provided  the  guess  is  good  enougn,  the  path  is  cor- 
rected by  making  the  integrated  value  of  H  smaller  and  smaller 
along  the  path. 
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In  a  given  problem,  the  degree  of  "goodness"  required  of  this 
nominal  path  when  using  the  Gradient  Method  corresponds  in  a  rough 
way  to  the  "goodness"  of  the  first  guess  for  U*(0)  when  using  extre- 
mals. In  the  former  case,  however,  one  can  be  guided  by  physical 
reasoning.  The  physical  meaning  of  the  value  of  U'(0)  is  by  no 
means  as  clear. 

In  the  preliminary  work  associated  with  this  thesis,  the  authors 
have  worked  on  a  means  to  generate  a  point  U'(0)  within  R  from  a 
nominal  path.  In  the  sample  problem  undertaken,  the  method  works 
well  and  provides  a  value  for  U'(0)  which  is  quite  close  to  U*. 
The  method  is  quite  simple  and  straightforward  and  adds  very  little 
complexity  to  the  computer  program  which  is  used  to  correct  from 
U»(0)  to  U*(0). 

2.3  The  fundamental  adjoint  set 

As  a  primary  tool,  we  make  use  of  the  fundamental  set  which 
is  used  by  Faulkner  ( 5  ).  The  fundamental  set  is  the  nxn  matrix, 


(2.1) 


U  = 


ull  u12 


.  u 


In 


U21  U22  •  '  '  U2n 


u  ,  u   .  .  .  u 
nl  n2       nn 


whose  value  at  t=0  is  the  nxn  identity  matrix  and  whose  elements 

are  chosen  so  that  u.  .  is  the  j'th  component  of  the  i'th  solution 

■•-J 
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to  the  adjoint  equation  (1.5).  Note  that  each  row  vector  of  the 
fundamental  set  is  a  particular  solution  to  the  adjoint  (1.5) • 
For  instance, 

(2.2)  ^(0)=  [^21u22.  .  .»g=  £100.  .  .0]. 

t=0 

How,  consider  a  particular  combination  of  the  row  vectors  U.  , 

(2.3)  U  =  cx  Ux  +  c2  U2  +  .  .  .  +  cn  Un, 

where  the  c's  are  constants. 
This  may  be  written  as 

(2.4)  U  =  CU,  where 

(2.5)  C  =  Jcx  c2  .  .  .  £}. 

Note  that  as  thus  defined,  U(0)  =  C,  so  that  the  C's  are  the  initial 
values  of  the  adjoint  variables.  At  any  time  t,  u.  is  the  ith  com- 
ponent of  U: 

(2.6)  \^\^l   +c2U2i+  '  •  '  +CnUni« 
or 

(2.7)  \  =  c%»  where 

u.  is  the  ith  column  of  the  fundamental  set. 

One  of  the  c's  may  be  chosen  arbitrarily.  If  some  state  variable 
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does  not  appear  explicitly  in  the  equations  of  motion,  then  the 
corresponding  adjoint  variable  is  constant „  If  this  state  vari- 
able is  to  be  maximized  then  it  will  be  convenient  to  make  the 
corresponding  c  equal  to  plus  one.  If  the  variable  is  free  at 
either  end  point  then  the  corresponding  adjoint  variable  is  zero 
everywhere,  and  a  different  c  will  be  chosen  arbitrarily.  In  any 
case,  we  must  be  sure  that  the  corresponding  adjoint  variable  is 
non-zero. 

ZA    Program  for  generating  the  C- vector  from  a  nominal  path 

In  this  section  a  program  is  given  for  generating  a  starting 
set  of  constants  U'.  A  nominal  control  program  is  chosen  from 
physical  reasoning  which  appears  likely  to  be  fairly  close  to 
optimal  and  which  will  generate  a  path  from  the  given  initial 
conditions  to  a  terminal  point  which  satisfies  or  nearly  satis- 
fies the  terminal  constraints.  This  may  be  simple  or  difficult, 
depending  on  the  problem.  If  such  a  program  cannot  be  found, 
then  this  method  for  starting  cannot  be  used,  nor  can  be  Gradient 
Method,  since  both  are  identical  up  to  this  point.  In  this  case, 
we  must  go  back  to  the  "guessing  game".  It  may  be  necessary  to 
"map"  the  U(0)  hyperplane  to  find  which  regions  give  good  start- 
ing points,  and  then  to  try  various  points  within  those  regions. 
That  is,  guess  the  c's  and  calculate  the  corresponding  trajec- 
tories in  some  systematic  way. 
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Suppose  however  that  a  satisfactory  nominal  control  program 
is  found.  The  equations  of  motion  are  now  integrated  numerically 
from  the  given  initial  conditions  using  the  nominal  control  pro- 
gram, stopping  when  a  terminal  constraint  is  met0  At  the  same 
time,  we  integrate  also  the  n  equations  of  the  fundamental  set. 
At  the  terminal  point  S(T),  each  of  the  elements  of  the  fundamental 
set  has  some  value  u. .(T). 

If  the  path  had  been  the  optimum  path,  equations  (1.1*0 
would  have  been  satisfied  at  the  point  S(T)6  We  will  choose  the 
vector  C  so  that  they  are  satisfied,  using  the  values  from  the  nom- 
inal path.  Then  this  vector  C  becomes  U'(0)  for  starting  the  cor- 
rective iterations.  If  the  program  does  not  converge,  the  nominal 
path  may  be  varied  somewhat  to  get  another  trial  C.  In  this  way, 
the  problem  of  guessing  the  initial  values  of  the  adjoint  variables 
may  be  eliminated,  as  in  the  Gradient  Method  and  still  the  optimum 
path  is  approached  through  extremals. 

It  is  to  be  noted  that  equations  (1.1*0  are  not  the  only 
relationships  which  are  used  to  furnish  the  required  number  of 
equations  in  the  C's.  According  to  the  particular  problem,  the 
Hamiltonian  may  be  constant,  or  a  particular  adjoint  variable  may 
be  constant.  These  allow  use  of  relations  calculated  for  t=0 
rather  than  at  t=T  as  above. 
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2.5  Faulkner's  method  of  differential  corrections „ 

A  method  was  given  above  for  obtaining  a  first  approximation 
for  the  starting  values  of  the  adjoint  variables  from  a  nominal 
path.  Let  us  now  suppose  that  such  a  set  is  in  hand  which  was 
obtained  in  this  way  or  in  some  other  way. 

If  the  equations  of  motion  (1.1)  together  with  the  adjoint 
equations  (1.5)  are  now  integrated  numerically  with  the  control 
chosen  according  to  the  maximum  principle  or  (1.7)  a  path  results 
which  is  an  extremal  in  that  it  furnishes  an  optimum  to  some  end 
state.  But  the  path  obtained  using  this  particular  set  of  initial 
conditions  does  not,  in  general,  satisfy  the  terminal  constraints 
(1.2)  and  the  transversal  conditions  (1.1^).  We  must  have  some 
means  for  correcting  these  constants  (which  we  have  variously  called 
the  C- vector  or  the  vector  U(0)  in  earlier  sections)  so  that  (1.2) 
and  (1.1^)  are  satisfied  at  some  point  S(T)„ 

The  method  used  here  is  the  method  of  differential  corrections 
developed  by  Faulkner  (ref.  5)  which  makes  use  of  differentials  in 
an  iterative  scheme  which  is  of  the  Newton-Raphson  type.  Making  use 
of  the  fundamental  set  and  the  development  of  Chapter  I,  the  method 
may  be  presented  quite  simply. 

Equation  (1.6)  provides  the  relation  (for  a  problem  without  cor- 

ners) 

T  r1 

JU&SJ      =      J      UF  <SPdto 
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For  every  state  variable  si  which  is  specified  at  t  =  0,  os^Co)  =  0, 
For  those  which  are  free,  we  will  require  Uj_(0)  =  0,  so  that  u.  ds^  =  0 
for  all  cases.  This  procedure  takes  care  of  separated  end  conditions 
in  a  simple  and  straightforward  manner „  This  case  causes  a  great  deal 
of  difficulty  with  some  methods  of  solution.  Note  that  there  are 
still  n  unknowns  at  the  start,  since  the  value  of  s^_  is  unknown  but 
the  value  of  the  corresponding  vu   is  known. 
With  this  simplification,  we  have 

m 

(2.8)  ftj  <5s]T  =  fo     (UFp)  6P  dt, 

where  the  dimensions  are  nxn  for  U„  nx1  for  oS„  nxm  for  F  and  mx1 
for  P.  From  (1.7)  we  will  choose  the  p's  so  that 

Hp  =  UFp  =  0. 
Now  if  S  and  t  are  fixed  and  the  c*s  are  varied  by  small  amounts  oCa 
then  P  must  vary  so  that 

(2.9)  6cUcFp  +  6PT  Hpp  =  0, 

where  the  superscript  T  indicates  the  transposed  matrix.  But  since 
U  =  CU  ,  then  Uc  =  U.  Making  this  substitutions,  and  noting  also 
that  H   is  the  symmetrical  mxm  matrix  whose  ij'th  element  is 

we  now  form  the  mx1  transposed  of  (2.9). 

(2.10)  (UFp)T  6cT  +  Hpp  6P  =  0 
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This  equation  is  solved  for  &P,  and  the  result  is  substituted 
into  (2.8). 


b?   =  -  (H  )"1(UF  )T  6cT 
PP     P 


(2.11) 


T  —1     T   T 

[j  fe]   =  -J  (UFp)  (Hpp)   (UFp)  c5c  dt 


T     o 
a  form  convenient  for  machine  programming. 

Now  suppose  each  term  in  the  integral  on  the  right  be  denoted 


as  I...  We  have  then  at  t  =  T 


(2.12) 


u 
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Since  one  c  is  arbitrary,  say  c-,  ,  then  6c-,   =  0  and  the  equation 
furnishes  n-1  relations  among  the  6c1  s  and  the  o's's. 

In  addition,  the  transversal  conditions  provide  n-k  equations 
whose  variational  form  is,  for  the  problem  to  be  considered  here, 

ou.  =  6cui4 
Now  for  corrections  we  let 


(2.13) 


6u±  =  -  u^T)  -  u±(T)  6t  =  6qu. 


and  also,  for  the  k  constrained  variables, 


6s±  =  -  s^T)  -  Si(T)  6t  +  sc£r) 
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where  s  (T)  is  the  given  terminal  constraint  for  the  ith  variable. 
Since  there  are  k  of  these  variables  constrained ,  there  remain 
(n-  k-l)  elements  &s^  in  equation  (2.12).  These  may  be  eliminated, 
yielding  k  equations.  The  n-k  equations  (2.13)  are  adjoined  to 
these,  for  a  total  of  n  equations  involving  the  unknowns  6t,  6c2, 

6Co,  .  .  .  6"Cn  in  terms  of  the  integral  elements  I^.„  the 
elements  of  U(T),  the  deriviatives  at  t=T  and  the  differences  (2.13) 
The  equations  are  solved  for  the  unknowns,  and  these  are  applied  as 
corrections  for  the  start  of  the  next  iteration. 
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Chapter  III 

DEVELOPMENT  OE  THE  PROBLEM 

3.1  General  description  of  the  problem  and  assumptions. 

The  problem  described  in  this  thesis  is  that  of  deter- 
mining the  optimal  angle  of  attack  program  p(t),  which  will 
maximize  the  distance  covered  over  the  earth's  surface  by  a 
hypersonic  glider  or  lifting  vehicle  which  has  been  injected 
into  some  initial  re-entry  point.   The  starting  point  of  the 
trajectory  is  specified  by  the  re-entry  injection  parameters, 
the  initial  conditions,  which  are: 

v{0)  =  25,920  feet/second  -  initial  velocity 

h(0)  =  30C ,000  feet       -  initial  altitude 
(3.1) 

x(0)  =  0  nautical  miles   -  initial  distance 

z(0)  =  0.18  degrees       -  initial  flight  path  angle 

These  parameters  were  taken  from  reference  (8)  so  that  a 

comparison  could  be  made  between  the  Gradient  Method  used  by 

earlier  authors  and  the  method  of  differential  corrections 

which  was  used  in  this  thesis.   Also,  the  value  of  wingload- 

ing,  the  value  of  acceleration  due  to  gravity,  the  value  for 

a  standard  earth  radius,  and  the  value  of  a  standard  nautical 

mile  were  taken  from  the  same  reference. 

Eingloading  =  m£0  =27=3  lb/ft2  ,  where: 

A 
m  =  mass  of  the  vehicle 

A  =  wing  plan  form  area 

g0  =  32.2  ft/sec  -  acceleration  due  to  gravity  at 

the  earth's  surface 

The  model  atmosphere  used  was  based  on  the  tables  of  reference 

(9).   ^he  atmospheric  density  was  assumed  to  have  the  form 

31 


p-oQ~      ,  where  the  parameter  b  was  calculated  to  fit  the  tab- 
ular values  at  h  =  0  feet  and  h  =  200,000  feet.   It  was  fur- 
ther assumed  that  the  atmosphere  was  spherically  symmetric  and 
fixed  with  respect  to  the  earth  for  simplifying  reasons, 
3.2  Equations  of  motion. 

The  coordinate  system  used  is  depicted  on  Fig.  3.1  and  the 
equations  of  motion  are: 

*  =  rTH  v  cos  z  ' 


(3.2) 


ft   =  v   sin  z  , 

D 
v  =  -  -  -  g  cos  z  , 

m   ° 

z  =  h  +  rth  cos  z  -  f  cos  z 

where:   R  =  3440  nautical  miles,  the  radius  of  the  earth,, 
D  =  |(/>v2A  CL)  , 

L  =  \{r2k   CD)  , 

r  R  -.2 
g  "  s0LH+hJ  ' 

The  coefficients  of  lift  and  drag  are  from  reference  (8)  and  • 
are : 

CT  =  Cr~  sin  d  cos  p  I  sin  p 

Li        JjU 

CD  .  CLL  lsin3pl+  CDQ 

Where:   CLn  =  1.82  ;   CDL  =  1.46  ;   CD0  =  0.04?  ;  and  p  is  the 
angle  of  attack  of  the  vehicle. 
3.3  Adjoint  equations. 

rihe  system  of  differential  equations  which  is  adjoint 
to  the  variational  equations  of  the  equations  of  motion  (3,2)  is 
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ux   =   0 


(3.3) 


r     Rv  i  rl   cD        cp  i 

uh  =    LrR7h72    cos   zjux  +    ^  ^  -  ^  sm  z]  u 

fv    cos    z        1_    cL   +  £g  cosz] 
UR+h)2      "my   eh        dh       v     Juz 

f-R  cos   zi  r    .        1  rl   6D7 

"v  =  I— HTB — K  "lsln  zjuh  +  li  evK 


mv 


2  -  75 


eL  l 


cv  mv 


COS     Z    _    g    COS    Z' 

R  +   h  v2        j " 


lh 


"  Rv   sin   zl  r  i  r  i 

uz  -  1  ~ETh— 'Jux  -  [v  cos  zJuh  +  [e  cos  ZJU, 


+    fv    rin  z    .   S   sin   z1u 

+    LR  +  n  v       Juj 


The  fundamental  set  of  adjoint  equations  U  or  equivalently 

(u,  .)  was  chosen  such  that  the  fundamental  set  at  time  t  =  0 
1 J 

is  the  identity  matrix  of  rank  four.   That  is: 

a-  <<v  = 


(3.^) 


h1 

X 

u2 

X 

"1 

4 

ux 

< 

< 

-2 

4 
uh 

u1 

V 

u2 

V 

«? 

4 
uv 

u1 

.  z 

2 
uz 

z 

uz 

i t   rT  ° 

and  further:   (u, ,) L  =70  (u,.)dt  ,  where  the  u. .  are  defined 
by  the  adjoint  equations  (3«3)«   We  have  the  further  relation 
that:  U  =  [u  u.  u  u  ]  =  [c-,  c2  c7  c^][u]   which  is  a  more 
convenient  form  for  computations  using  Fortran  programming. 

Tn  general,  the  solution  of  the  adjoint  equations  is 
only  determined  to  within  a  multiplicative    constant  ,  hence 
we  can  choose  one  of  the  c's.   Tn  this  problem,  we  chose  c,=l„ 

Tf  we  consider  the  equations  of  motion  in  vector  or 
matrix  form,  we  have  the  following  relationships: 
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s 


x 

h 

v 

.z 


F  =  S  ; 


0 
0 

m 

mv 


-*    c  ^ 

where  P„  =  ■3— 

P    cp 


The  Hamiltonian  is  given  by  H  =  U°F. 

3.^  Maximum  principle. 

For  any  extremal,  we  must  satisfy  the  condition  that 

U»F  =  0  for  any  point  on  the  extremal.   This  condition 

uniquely  determines  the  angle  of  attack  program  for  that 

particular  extremal  since: 

V'%  =   §|u  -v~u  =  0 
P   cp  z    dp  v 

Substituting  for  •*—  and  -g—  and  solving  for  p  we  obtain? 


(3.5) 


p  =  arctan 


3CDL  v  1 
L2CLQ  uzJ 


l\l 


8 


CL0 

C~~uz 

LDL  Z 


"21 


-  u. 


F 


f  or  -  j  <   p  <  j 

?.5   Fnd  conditions. 

Since  we  wish  to  maximize  the  total  distance  traveled 
over  the  earth's  surf see,  we  must  have  the  following  conditions 
at  the  end  time  t  =  T: 

x(T)  =  maximum 

(3.6)  h(T)  =  0 

v(T)  ,   z(T)  ,   and  T  free. 
3.6  Variational  equations. 

On  any  curve  without  corners,  we  have  the  following 
relation  between  the  adjoint  and  the  variations: 

T 

(3.7)  [u^6x  +  u^&h  +  u^cv  +  u^6z]T=^  "u'fe  cp  dt   j  =1,2,3,4 
For  an  extremal  we  must  also  satisfy  the  Euler  equations  or 
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the  maximum  principle  of  Pontryagin  so  that  we  have 

rfi'_.  _*  _»  _* 

J  ~U*F  dt  =  maximum  and  U'Fn=  0  along  the  extremal.   If  we 

consider  x,h,v,z,  and  t  fixed  and  change  the  constants  then 
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ec\6c2  +  dc„6c3  +  fo  ec^J,pp  +  u'p™ &p 


3 


PP 


-g~  =  U   ,  1  =2,3,^,  we  can  solve  for  6p  and  substitute  in 
equation  (?.7)  which  will  give  us     four  equations  with 
which  we  can  correct  the  initial  values  of  the  c's  for  a 
terminal  variation  of  the  altitude  h: 


(3.8) 


:,,-! 


u^&x  +  u£&h  +  u^6v  +  u^zj  r   -51  /  (tf^p)(U-Fp)dt6c 
<v     —^  .  _»1_,  .  PP 


Let  I.  ,= 

J  "0 


(UVF   (U-P  ) 

■ *- H-dtec. 


U.P 


PP 

2     3     2i 
Since  u  =  1 ,  u  =  u-  =  uv  =  0  for  all  t,  and  noting  that 

X  XXX 

I. .  =  I**,  we  can  solve  the  following  set  of  equations  for  fch 
in  terms  of  the  6c, : 


'4 

2 

uv 

21 
uz 

*6h" 

T     T    T 

22   23  2k 

■2 

< 

a3 

z 

ov 
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23  33  3*4- 

k 

k 

k 

uh 

uv 

uz. 

T 

6z 

.I2^-  I3^  T^ 

6c2 

fc. 

T 

.6c4. 

(3.9) 


Prom  the  transxersal  conditions,  we  have  that; 
uy(T)  =  0 
(3.10)   uz(T)  =  0 

(xuy+huh)T  =  0 
^rom  the  last  equation  of  (3.10)  we  have:  uh(T)  =  (-cot  z).T 

-A   _> 

Since  the  Hamiltonian  H  =  U«F  and  H(T)  =0  from  equations 

(3.10),  then  the  Hamiltonian  is  identically  zero  for  all  t. 

This  fact  will  be  used  to  determine  one  of  the  constants  (c2) 

and  effectively  reduce  the  problem  of  guessing  the  initial 

c's  to  start  the  calculations. 

35 


The  transversal  conditions   (3.10)  give  us  three  more 
independent  variational  equations  associated  with  the  end 
coni  it ions : 

6uv(T)  =  [u^6c2  +  u^6o3  +  v^bck]T 

(3.1D   6uz(T)  =  [u^6c2  +  \xhc3   +  u^6o4]T 

6[xux+  huhlT  =  h(T)[uh6c2  +  u^6c3  +  uh6cJ!f]T 

3.7  Differential  corrections 

The  method  of  differential  corrections  given  by 
Faulkner  in  reference  (5)  was  used  for  the  solution  of  the 
optimization  problem.   In  general,  we  will  guess  the  initial 
values  associated  with  the  adjoint  equations  and  the  terminal 
time  T.  For  this  problem,  we  must  choose  03,  c-a ,  c^ ,  and  T. 
We  then  correct  these  parameters  by  calculating  an  extremal 
with  the  initial  guess  and  then  solve  for  corrections  to  the 
parameters  in  terms  of  the  end  values  associated  with  the 
extremal.   For  all  end  conditions  of  the  form  s(T)  =  S ,  we  will 
generally  have  some  error  in  s(T).   We  set  6s (T)  =  S  -  s(T)  - 
s(T)  6T.   For  this  problem,  we  set  ch(T)  =  h(T)-h(T)£T  and  in 
a  similar  manner,  we  obtain  from  equation  (3.11)° 

-uv(T)  =  [u^6c2  +  u^6c3  +  u^£cJT  +  uv(T)cT 
(3.12)   -uz(T)  =  [uz6c2  +  u^£c3  +  Ac^]T  +  uz(T)fT 
[cot  z]T  =  h(T)[u^6c2  +  u^ec3  +  uh£c,JT 

+[xux  +  huh  +  huhlT6T 

These  equations  together  with  the  solution  to  (3. 9)  for  6h  in 
terms  of  the  6c' s  will  give  a  set  of  equations  from  which  we 
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solve  for  the  fc's  and  fT  to  make  corrections  In  the  starting 
values  for  the  subsequent  trajectory. 

3.8  Computation  of  approximate  starting  values  by  use  of 
nominal  trajectories. 

In  general,  for  a  nth  order  set  of  differential 
equations,  one  must  choose  (n-1)  of  the  constants  of  the  sol- 
ution to  the  -sdjoint  equations  to  start  the  oroblem.   This  is 
generally  a  difficult  problem  in  itself,  since  the  constants 
are  not  related  to  the  ohysical  aspects  of  the  problem  in  any 
discernable  manner.   One  of  the  primary  aspects  of  this  thesis 
was  to  find  a  method  which  could  be  used  to  find  a  logical 
choice  for  the  initial  values  of  these  constants.   Fortunately 
the  Hamiltonian  was  identically  zero  in  this  problemc   This 
fact,  together  with  the  transversal  conditions  ( 3»1°)  enabled 
us  to  solve  for  approximate  constants  by  first  computing  a 
nominal  trajectory;  here  we  choose  a  nominal  angle  of  attack 
program  p(t),  which  seemslikely  in  a  physical  sense  to 
give  a  maximum  distance  trajectory.   The  first  such  nominal  angle 
of  attack  program  tried  was  to  use  the  maximum  lift  over  .irag 
ratio  where  p=20.5  degrees,  a  constant.   This  did  not  work, 
We  then  used  the  following  angle  of  attack  program  for  which 

we  are  indebted  to  Professor  Faulkner" 

.0 


P(t)  =  5k. 7 


1  -  h(t) 

750,000  J 


The  value  of  5^.7*  for  p  is  the  value  at  which  we  have 
maxium  lift.   It  was  supposed  that  the  ar.gle  of  attack  was 
never  negative  and  from  the  data  of  constant  angle  of  attack 
programs,  the  value  of  75O,OO0feet  was  chosen  to  insure  that 
p  was  never  negative. 
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solve  for  the  fc's  and  £T  to  make  corrections  in  the  starting 
values  for  the  subsequent  trajectory. 

3.8  Computation  of  approximate  starting  values  by  use  of 
nominal  trajectories. 

In  general,  for  a  nth  order  set  of  differential 
equations,  one  must  choose  (n-1)  of  the  constants  of  the  sol- 
ution to  the  adjoint  equations  to  start  the  oroblem.   This  is 
generally  a  difficult  problem  in  itself,  since  the  constants 
are  not  related  to  the  ohysical  aspects  of  the  problem  in  any 
discernable  manner.   One  of  the  primary  aspects  of  this  thesis 
was  to  find  a  method  which  could  be  used  to  find  a  logical 
choice  for  the  initial  valves  of  these  constants,   Fortunately 
the  Hamiltonian  was  identically  zero  in  this  problem.   This 
fact,  together  with  the  transversal  conditions  ( 3«10)  enabled 
us  to  solve  for  approximate  constants  by  first  computing  a 
nominal  trajectory;  here  we  choose  a  nominal  angle  of  attack 
program  p(t),  which  seems  likely  in  a  physical  sense  to 
give  a  maximum  distance  trajectory.   The  first  such  nominal  angle 
of  attack  program  tried  was  to  use  the  maximum  lift  over  drag 
ratio  where  p=20.5  degrees,  a  constant..   This  did  not  work. 
We  then  used  the  following  angle  of  attack  program  for  which 

we  are  indebted  to  Professor  Faulkner : 

.0 


P(t)  =  5^.7 


1  -  hill   __ 
750,000 


The  value  of  5^*7     for  p  is  the  value  at  which  we  have 
maxium  lift.   It  was  supposed  that  the  ar.gle  of  attack  was 
never  negative  and  from  the  data  of  constant  angle  of  attack 
programs,  the  value  of  750, 000 feet  was  chosen  to  insure  that 
p  was  never  negative. 
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With  this  choice  for  the  angle  of  attack  program,  we 
obtained  a  good  set  of  starting  constants  by  integrating  the 
fundamental  set  U  along  with  the  equations  of  motion  until  h  =  0 
At  this  time,  we  solve  the  set  of  equations: 


(3.13) 


for  C2 ,  Co,  and  c/j,.   These  equations  are  the  transversal 

equations  (3 .10) .  This  problem  of  guessing  the  cBs  is  further 

reduced  due  to  the  Hamiltonian  being  identically  zero,  since 

for  a  given  choice  of  Co  and  c^ ,  the  angle  of  attack  is  uniquely 

determined  from  equation  (3*5)  and  then  we  may  solve  for  C2  from 

the  Hamiltonian  at  time  t  =  0.   If  we  have  the  case  where  h(0)= 

0,  then  we  have  a  singularity  and  connot  solve  for  Co ,  but  in 

this  case,  the  choice  of  03  is  arbitrary.   Effectively,  this 

problem  is  then  reduced  by  an  order  of  one  and  we  need  only 

guess  two  constants,   We  felt  that  if  we  had  an  extremal,  we 

would  determine  its  constants  in  this  manner  and  hence  we  could 

which 
get  estimates  for  the  constants  from  a  trajectory^ was  a  "good" 

approximation  to  an  extremal. 

3.9  Computational  scheme. 

This  problem  was  orogrammed  on  a  CDC  1604  computer 

using  a  Runge-Kutta  integration  method.   The  computational 

scheme  is  as  follows: 

i)   We  need  a  nominal  trajectory  to  start.   To  get  it,  we 

first  intergrate  the  equations  of  motion  0*2)  and  the  fundamental 

set  (3.4)  using  a  programmed  angle  of  attack  and  ending  the 
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trajectory  when  h  =  0.   At  this  time,  we  solve  for  the 
initial  constants  c2 ,  Cj ,  and  c^  from  equation  (3,l}\     We  then 
"se  these  constants  to  start  the  first  optimal  trajectory, 

ii)   Then  we  compute  the  first  optimal  trajectory  where 
the  angle  of  attack  program  is  determined  by  the  maximum 
principle.   We  integrate  the  equations  of  motion,  the  fundament- 
al set,   and  the  nine  I1j  in  equation  (3.9)  a  total  of  twenty- 
five  equations.   We  terminated  the  integration  when  one  of 
the  following  conditions  is  met:   h  =  0,  u  =  0 ,  u.=  0 ,  or 

I   !  ^  TT 

lzl=  p»  since  it  was  felt  that  the  end  of  the  significant  part 
of  the  trajectory  had  been  reached.   We  felt  that  in  this 
problem  no  optimum  would  occur  for  u  <  0  since  the  resulting 
high  drag  would  dissioate  energy  needlessly.   The  supposition 
that  this  is  a  maximum  distance  trajectory  dictates  that  I z I <   ~« 

Since  we  do  not  have  t  appearing  explicitly  in  the  equat- 
ions and  since  the  terminal  conditions  determine  T,  we  do  not 
use  6T  in  our  calculations.   We  only  need  the  6c' s  the  correct 
the  c's  to  start  the  next  trajectory. 

iii)   We  then  repeat  the  computations  as  in  the  above 
paragraph,  using  the  new  c's  for  the  next  trajectory.   We 
test  for  convergence  at  the  end  of  each  trajectory: 
(3.U)   If  [h(T)j2  +  [xux  +  huh]2  +  [uv(T)]2  +  Luz(T)]2  =£ 

we  say  that  we  have  converged  to  the  solution. 
3.10  Results  and  conclusions. 

The  optimum  trajectory  that  was  calculated  is  deDlcted 
on  Fig. 3.3  and  the  corresponding  angle  of  attack  program  is 
depicted  on  Fig.  3A  On  these  two  figures,  a  comparison  has 
been  made  between  the  results  of  this  thesis  and  results  from 
reference  (8)  ,  which  were  obtained  by  the  Gradient  (or  Steep- 
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est  Ascent)  Method.   The  comparison  can  r.ot  be  considered  as 
an  exact  one  since  the  integration  method,  the  time  step 
ofc  integration  interval,  and  the  approximation  to  the 
atmospheric  models  may  not  be  identical.   However,  there 
appears  to  be  a  significant  increase  in  the  maximum  distance 
obtained  using  the  method  of  differential  corrections. 

The  convergence  criteron  (3.  14)  we  initially  chose  turned 
out  to  be  a  poor  choice.   We  could  not  determine  how  close  we 
were  to  the  optimal  solution.   We  then  proceeded  to  "map" 
the  plane  of  Co  and  c/4,  by  using  a  systematic  choice  of  03  and 

Cj,  and  computing  the  corresponding  extremals <>   The  results  of 
the  mapping  procedure  are  given  on  Fig. 3.5,  Fig. 3.',  and  Figo3c7<> 
Fig. 3.5  shows  the  contours  of  distances  obtained  by  extremals 
as  a  function  of  c^  and  Cj,.   Fig.  3*6  gives  detailed  contours 
around  the  optimal  solution.   Pig. 3,7  is  a  cross  section  of  the 
"ridge"  which  occurs  in  the  c^  and  Cj,  plane.   It  is  a  con- 
jecture of  the  authors  that  some  difficulties  are  encountered 
by  a  gradient  method  solution  when  there  exists  such  a  ridge 
in  a  problem,  since  the  slope  is  nearly  zero. 


The  major  problem  encountered  was  the  apparent  instability 
of  the  problem  near  the  end  of  an  extremal,   We  attempted  to 
integrate  backwards  from  the  terminal  conditons  and  found  that 
this  was  virtually  impossible.   The  problem  seems  to  be  unstable: 
we  could  not  even  integrate  backwards  with  constant  p.   Back- 
ward integration  was  very  good  if  we  started  backward  from  a 
point  that  was  not  too  close  to  the  end  of  the  trajectory. 

Another  main  problem  encountered  was  that  of  determining 


the  end  time  on  a  trajectory.   We  first  thought  that  all  we 
need  consider  is  the  time  when  h  =  0  .   However,  we  found 
thai  other  conditions  effectively  terminated  the  useful  part 
of  the  extremal  earlier.   For  example,  most  extremals  led  to 
a  point  where  izi>  •=•• 

Theoretically,  on  the  maximum  distance  trajectory,  all 
terminal  conditions  would  be  met  simultaneously.   Tn  the  actual 
computations,  this  did  not  occur.   The  most  probable  reason 
that  we  have  this  apparent  discrepancy  is  that  the  round  off 
errors  in  the  calculations  and  the  inherent  error  in  the 
integration  method  are  large  enough  to  nrevent  us  from  satis- 
fying all  conditions  at  the  end  time, 

Tn  the  final  analysis,  we  found  that  a  program  which 
reduced  the  magnitude  of  the  corrections  would  lead  to  con- 
vergence.  The  method  used  was  to  calculate  a  set  of  c's 
from  a  nominal  trajectory  and  then  calculate  the  first  tra- 
jectory and  the  6c 's  corresoonding  to  that  trajectory , 
store  x,  the  c's,  and  cc's  of  this  trajectory  and  then  cal- 
culate the  next  trajectory.   If  the  distance  of  this  trajectory 
is  greater  than  the  distance  of  the  preceding  one,  we  proceed 
with  the  iterative  process.   Tf  the  distance  is  less,  then 
we  reduce  the  corrections  by  one  half  and  calculate  the  tra- 
jectory again.   If  we  still  have  less  distance,  we  further 
reduce  the  corrections  until  we  obtain  a  greater  distance. 
The  reason  we  developed  this  method  was  because  the  magnitudes 
of  the  fc's  were  of  the  order  of  magnitudes  of  the  c's 
was  noted  that  the  integrals  relating  6h  to  the  £c's  are 
improper  integrals.   This,  in  addition  to  the  unstable  char- 


acter  of  the  end  of  the  trajectory,  tended  to  make  the 
magnitude  of  the  6c*s  large;  however,  the  program  did 
correct  in  the  right  direction.   These   integrals  seem  to  be 
divergent  for  the  extremals  which  satisfy  the  transversal  con- 
ditions.  A  formulation  of  the  oroblem  suggested  by  Professor 
Faulkner  which  would  avoid  these  divergent  Integrals  by  using 
a  set  of  u's  which  were  determined  at  the  end  time  was  pro- 
grammed.  However,  this  program  was  not  checked  out  or  run 
since  the  original  integrals  did  not  actually  diverge  in  cal- 
culations made  on  the  computer.   The  method  suggested  was  to 

*  -it- 

use  a  fundamental  set  U   such  that  U  (T)  =  6,  .    After  com- 
puting an  extremal  using  U  ,  we  solve  for  the  matrix  of  con- 
stants ra,  .]  which  satifies  the  relation:   U(T)  =  [a44]U(T)« 
i  J  —        l  j  — 

We  then  compute  the  same  extremal  using  U  =  U,  +  cU~  ,  where 

-»    4-     -i 

U4  =  2_  a-  <u  »  i=  1»2.   We  then  have  the  relation 
.1  =  1   1J 

2   '  v>'   dt  to     and  we  solve  for  6c  to  correct  the 
-_*-_ — 

u   U  -F 

PP 

starting  va'jues.   Further,  this  integral  is  not  divergent. 
The  first  integration  method  used  in  calculations  was 
rectangular  integration  and  it  was  found  that  trajectories 
thus  calculated  were  very  rough.   A  comparison  of  the  rectang- 
ular integration  method  and  the  Runge-Kutta  method  was  made 
and  we  chose  the  latter  method  due  to  its  accuracy  and 
smoothness. 
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